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We present a systematic overview of granular deposits composed of ellipsoidal particles with different parti- 
cle shapes and size polydispersities. We study the density and anisotropy of such deposits as functions of size 
polydispersity and two shape parameters that fully describe the shape of a general ellipsoid. Our results show 
that, while shape influences significantly the macroscopic properties of the deposits, polydispersity plays appar- 
ently a secondary role. The density attains a maximum for a particular family of non-symmetrical ellipsoids, 
larger than the density observed for prolate or oblate ellipsoids. As for anisotropy measures, the contact forces 
show are increasingly preferred along the vertical direction as the shape of the particles deviates for a sphere. 
The deposits are constructed by means of an efficient molecular dynamics method, where the contact forces are 
efficiently and accurately computed. The main results are discussed in the light of applications for porous media 
models and sedimentation processes. 
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I. INTRODUCTION 

The granular character of many processes in nature has mo- 
tivated the study of granular materials and in particular granu- 
lar packings since Kepler's time[l |. One of such processes is 
the deposition of particles subject to gravity 12, which under- 
lies the formation of sandstones, ceramics and in general, the 
emergence of porous materials. Modelling porous media||4][5) 
in a realistic way is important for instance to understand per- 
meability of soils for petroleum engineering or in material sci- 
ences. 

However, a model for grain deposition should incorporate 
two ingredients not always easy to combine. First, empiri- 
cal studies [6 1 show that the shape of grains significantly de- 
viates from the sphere with a non-negligible polydispersity in 
their sizes. Until recently most computational studies dealt 
with round particles for simplicity, in solving specific prob- 
lems such as obtaining the optimal packing!?] 0> construct- 
ing space-filling configurationsll9l 4TTI or studying elongated 
shapes[ 3]. Some progress was recently made when addressing 
the influence of particle shape on the macroscopic observables 
of the agglomerate, focusing on ellipsoidal particles lfT2l [T3l 
and even in more general shapes lfT4l[T5ll , without accounting 
for the size polydispersity. 

Second, deposition processes are fundamentally sequential 
and therefore should be modelled through sequential proce- 
dures. The procedure introduced by Donev et a/|Q] [12] for 
jammed configurations, captures some of the general features 
observed in granular deposits, such as their density and av- 
erage coordination number, although it is based on a non- 
sequential procedure. However, deposition under gravity does 
not generally lead to a jammed state. Recently, some sequen- 
tial procedures were already carried out in two-dimensional 



deposits of mono-disperse elongated particles [2] E). 

In this paper we focus on three-dimensional polydis- 
perse deposits of general ellipsoids and study how macro- 
scopic properties of such deposits depend on the shape and 
size polydispersity of particles. We show that the density 
and coordination number of deposits behaves similarly to 
what is observed in jammed systems. However, in contrast 
with randomly constructed jammed systems, deposits become 
strongly anisotropic when the shape of the particles deviate 
from the sphere. Furthermore, we find that the size polydis- 
persity plays a minor role when compared to the shape. 

To this end, we develop an efficient algorithm for the calcu- 
lation of the interactions forces at all contact points within 
a three-dimensional deposits of polydisperse ellipsoids and 
use it in a Molecular Dvnamics |fT6l WT\ method to perform 
our simulations. Further, we introduce two shape parameters 
which describe any ellipsoid, enabling a systematic study of 
the role of particle shape in granular deposits. 

The outline of the paper is as follows. In Sec.[ll]we describe 
the algorithm for calculating the contact forces between col- 
liding ellipsoids with arbitrary size and shape as well as the 



simulation setup. In Sec. Ill the control parameters describ- 



ing the polydispersity in the size and the shape of a general 
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ellipsoid are introduced. In Sec. IV the dependency of macro- 
scopic observables taken from the entire deposit are studied as 
functions of the control parameters. We focus on the density, 
coordination number and anisotropy of the deposits. Discus- 
sion and conclusions are given in Sec.[V] 



II. MODELING THE CONTACT FORCES BETWEEN 
ELLIPSOIDAL PARTICLES 

When dealing with large systems, calculating the inter- 
particle forces is by far the most time consuming part of the 
computations. This is due to the fact that, in general, every 
particle can interact with every other one. In granular materi- 
als, however, the inter-particle forces are short range making 
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FIG. 1: (Color online) Sketch of the contact vector £ used to define 
the normal repulsive force between two overlapping ellipsoids A and 
B. The contact point is chosed to be X c after properly deriving points 
and Xm at me surfaces of the ellipsoids (see text). The full 
procedure is sketched in Fig. [2] 



techniques such as verlet list and linked-cell algorithm indis- 
pensable in reducing the computation time. These technique 
are employed to eliminate the redundant calculations for pairs 
which are too far apart to have any interaction |[T6ll . 

In time-driven Molecular Dynamics (MD) methods of gran- 
ular materials the system is evolved for one time step allowing 
the particles to overlap. Then normal repulsion forces on the 
particles are calculated as a function of the amount of over- 
lapping, which is defined usually by a line segment. Here we 
refer to this line segment as the contact vector (see Fig. [TJ. 

The contact vector between two ellipsoids can be defined in 
more than one way. The common method, introduced by Per- 
ram and co-workers 1 1 8 1 and developed by others |[T3~1 [T9l l20l . 
is based on finding two points, one on the surface of each ellip- 
soid by minimizing a potential through a variational problem 
with one geometric constraint. The two obtained points can be 
used to define the contact vector. However, the contact vector 
defined this way is not necessarily normal to the surface of the 
particles, particularly when particles are strongly asymmetri- 
cal, leading to non-physical motions of the particles. 

An alternative method consists in finding two points on the 
surfaces of the ellipsoids at which the corresponding surface 
normals point exactly in opposite directions. Consequently, 
the vector connecting these two points is normal to the sur- 
faces of both ellipsoids and can be used as a well-defined 
contact vector. This method is, however, much more com- 
putationally expensive than the former one. See Ref. 12D for 
more details. 

To retain the computational efficiency, we use the former 
method and introduce an additional step for correcting the di- 
rection of the contact vector. 

The geometric potential of the ellipsoid A is defined as: 

f A (X) = (X- X A ) T A(X - X A ), (1) 

where X and X A are the Cartesian coordinates of a point 
on the surface and the centroid of the ellipsoid, respectively. 
The equation f A (X) — 1 describes the surface of ellipsoid 
A. Considering a pair of ellipsoids A and B, two minimum 
points xffl and Xm are found on the surfaces of A and B by 



minimizing f B and f A , respectively. These points are shown 
schematically in Fig{T] To obtain the minima of f A (X) sub- 
ject to the constrain of being on the surface of ellipsoid B we 
minimize the following auxiliary potential: 

f(X) = (X-X A ) T A(X-X A )+\ [(X ~ X B ) T B(X - X B ) - 1] , 

(2) 

where A is the Lagrange multiplier associated with the con- 
strain. Setting to zero the gradient of this function with respect 

toJT, 

V/pO - A(X m - X A ) + XB{X m - X B ) = 0, (3) 

we obtain an equation for optimum points X m for A: 

X m = (A + AB)- 1 (AX^ + XBX B ). (4) 

To derive an equation for A we left-multiply Eq. ^ by 
(X m — X B ) T and use the fact that X m lies on B and that the 
gradients of the potential functions of the ellipsoids at mini- 
mum point are in opposite directions: 

A = | (X m - X B ) T A(X m -X A )\. (5) 
Equations (Wi and (Dj) can be solved iteratively to find Xm up 

— (B) 

to the desired precision. The second point X m can be found 
similarly. 

The contact point X c where a contact force is exerted on 

( A \ 

both ellipsoids is defined as the midpoint between X m and 
Xm ^ ■ The contact vector for ellipsoid A is defined asj^i = 
X { m ] - X ( ® ] and for ellipsoid B, £ B = -£4 (see Fig.[l}. 

Although Eq. Q involves inversion of a 3 x 3-matrix, find- 
ing the minima can be very efficient by choosing the initial 
points properly. Good approximations for minima are the 
minima from the last time step or, if the ellipsoids were dis- 
joint in the last time step, the closest points on the surfaces of 
the ellipsoids which are known from collision detection pro- 
cedure described later in this section. These points are then 
used as initial points, reducing the number of iterations sig- 
nificantly. 

The overlapping vector calculated in this way, however, is 
not necessarily normal to the ellipsoids at the contact point. 
In fact, the normal directions to A and B at points X&P and 

(B) 

X m deviate significantly from each other when the particles 
are strongly aspherical. To overcome this problem we intro- 
duce an additional step for correcting the contact vector, as 
described next. 

If the calculated contact vectors deviate more than a given 
amount from the gradients of the ellipsoid's potential at the 
contact point X c , a new direction is calculated by averaging 
the gradients at X c . Then, two new points are calculated from 
the intersection of the line along this direction passing through 
the contact point. These new points are used to define a new 
contact vector. Although it is not generally possible to define 
a vector which is along the gradients of both ellipsoids at the 
contact point, this approximation gives satisfactory results in 
terms of energy conservation. 

Having described the procedure for calculating the contact 
of two colliding ellipsoids, we next need an efficient way to 



3 




Compute 
the minimum points, 
the overlapping vector 
and 

the contact point 



FIG. 2: (Color online) Flow-chart of the algorithm for detection and 
calculation of the contact points between two arbitrary ellipsoids. 



determine when two previously disjoint ellipsoids collide. To 
this end, we use the technique proposed by Wang et q/ l25l 
which is described briefly as follows. 

In homogeneous coordinates the equation for the surface of 
a general ellipsoid is given by 



A : X T AX = 0, 



(6) 



where X is any point on the surface of the ellipsoid and A is 
the 4 x 4-matrix of the affine transformation which transforms 
the unit sphere centered at the origin into the ellipsoid. 

Given two ellipsoids A : X T AX = and B : X T BX = 
0, their characteristic equation is defined as f{y) — det(^A + 
B), whose roots are the eigenvalues of matrix — A _1 B. Fur- 
thermore, it can be shown 11251 that (i) f(y) — has at 
least two negative roots, (ii) f(v) = has two distinct pos- 
itive roots if and only if the two ellipsoids are disjoint, (iii) 
f(v) = has a positive double root if and only if the two el- 
lipsoids are externally touching. Therefore, by examining the 
eigenvalues of A _1 B it can be determined if two ellipsoids 
collide. 

When the ellipsoids are disjoint, the corresponding eigen- 
vectors give the coordinates of the four vertices Vi, i = 



1, . . . , 4, of a tetrahedron which is self-polar for both ellip- 
soids A and Bll26l i.e. V^AVj = V?BVj = 0,Vi ^ j. 
Two of the vertices lie outside both ellipsoids (external ver- 
tices), while the other two are contained inside the ellipsoids, 
one in each (internal vertices). In this case, a set of separation 
planes can be defined using three points exterior to both ellip- 
soids, two of them being the external vertices and the other a 
point on the line segment connecting the internal vertices. We 
choose the third point as the midpoint between the surfaces of 
two ellipsoids on this line segment. 

Having the separation plane between the two disjoint ellip- 
soids, in the subsequent time step we check if the ellipsoids in- 
tersect that plane. If not, the ellipsoids are still disjoint and the 
procedure stops. This helps to eliminate the need for checking 
for collision between the ellipsoids which is more computa- 
tionally expensive. Only when any of the ellipsoids intersects 
the separation plane it is necessary to check again for the colli- 
sion between the ellipsoids. If the ellipsoids turn out to be still 
disjoint their separation plane and closest points are updated 
to be used in the next time step. If the ellipsoids do collide 
the contact vector is calculated and saved to be used for the 
calculation of the contact force. 

Figure [2] shows the flow-chart of the complete procedure 
which is performed for all potentially overlapping pairs of el- 
lipsoids. It is worth noting that, to further speed up the pro- 
cedure, only the pairs whose spherical envelopes intersect are 
evaluated. 

The contact force can be decomposed to normal and tangen- 
tial forces. The tangential forces such as static and dynamical 
frictions are generally derived from the normal force. For cal- 
culating the magnitude of the normal force we use Hertzian 
modellfl6ll: 



F n CX e 



3/2 



dt' 



(7) 



with £ = |£| being the length of the contact vector between 
the particles and K being the dissipative constant depending 
on the material viscosity [22]. Since we are only concerned 
with the properties of the deposit at rest, we choose K suffi- 
ciently large for the deposit to relax rapidly. This relaxation is 
enhanced by including a tangential force F t — /iF n v t , fi be- 
ing the dynamic friction coefficient and v t the relative tangent 
velocity of the particles. 

Finally, the procedure described above is used in a MD sim- 
ulation together with standard methods, namely a prediction- 
correction integrator for solving the equations of motions, 
quaternions for describing the orientation of the particles and 
linked-cell algorithm to eliminate unnecessary collision de- 
tections. For details of the MD methods see e.g. Ref. |[T6l and 
references therein. 



m. CONTROL PARAMETERS FOR SIZE AND SHAPE 

The size and shape of the particles are the main parameters 
whose effects on the properties of the system are to be studied. 
In the following we explain how these are defined and chosen 
when preparing the samples. 



FIG. 3: Volume distributions V(r) for different values of a and j3 as 
defined in Eq. ([8}. The inset shows the corresponding size distribu- 
tions n(r) (see text). 



We define the size of an ellipsoidal particle as the radius of 

1 

the sphere with the same volume, i.e. r — (abc) 3 with a, b 
and c being the three semi-axis radii of the ellipsoid. We refer 
to a sample as mono-sized or monodisperse in size when all 
its constituting particles have the same volume, i.e. the same 
size r. To introduce size polydispersity we adopt the approach 
by Voivret et al &3\ : instead of directly choosing the size dis- 
tribution n(r) of the particles, we consider the distribution of 
the total volume V(r) — ^irn(r)r 3 of all particles with size 



We chose V(r) as 



V(r) = 



C 



(r—r m i n ) c 





if r„ 



< r < r n 
otherwise, 



(8) 



where a and j3 determine the shape of the distribution and C 
is the normalization factor. Figure [3] shows V(r) for differ- 
ent values of a and f3 and the corresponding size distributions 
n(r) (inset). For a = f3 this function takes a symmetric shape 
with a peak at (r m i n + r max )/2. For the samples studied in 
this work we have chosen a = j3 = 3 which corresponds to 
a distribution (solid line in Fig. [3} very close to a truncated 
Gaussian. 

The width of the distribution is controlled by the polydis- 
persity parameter <5 r which is defined as 



8 r = 



(9) 



For S r = the particles are monodisperse while 6=1 corre- 
sponds to infinite polydispersity. 

The shape of an ellipsoid is characterized by two parame- 
ters here defined as r\ = a/b > 1 and £ = b/c > 1, with 
a > b > c. For prolate (a = b < c) and oblate (a = b > c) 
ellipsoids, the shape can be fully characterized by the aspect 



FIG. 4: (Color online) A deposit consisting of TV ~ 4200 ellipsoids 
with shape parameters r\ — 1.4 and £ = 1.8 a polydispersity of 
S r = 0.6, that is, r min = 0.016 and r max = 0.064 (see text). 
Colors are arbitrarily chosen for better visualization. 



ratio a/c = 7/£. The particular case r\ = £ = 1 corresponds 
to the sphere. These shape parameters together with the size 
r > fully specify any ellipsoid. 

Similar to the size polydispersity, two additional parameters 
5 n and 5q can be similarly defined for the polydispersities in 
the shape of the particles. In this study we consider systems 
of mono-shaped particles, describing each deposit by 77, ( and 
size polydispersity 5 r . 

In this work, all the deposits are generated by releasing par- 
ticles with randomly chosen positions and orientations and 
letting them fall in an open box, under gravity along the ver- 
tical direction. The box is limited from below by a rigid wall 
(ground) and is periodic in x and y directions. In order to 
reduce the boundary effects the first few bottom layers of par- 
ticles are not considered in our analysis. Figure [4] shows a 
deposit of ellipsoidal particles with 77 = 1.4, ( = 1.8 and 
S r = 0.6. 



IV. MACROSCOPIC PROPERTIES OF THE DEPOSIT 

The macroscopic properties considered in our study are 
the packing density, the average coordination number and the 
anisotropy. 



A. Packing density and coordination number 

Figure [5] shows the packing density in deposits of N ~ 
3000 mono-disperse ellipsoids as a function of the shape. The 
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FIG. 5: (Color online) (a) Density p of monodisperse deposits as a 
function of the shape parameters r\ and £. The density is approxi- 
mately symmetric with respect to the diagonal 77 = £ where it attains 
a maximum, (b) The density as a function of shape along the diag- 
onal T) = £ (solid line) together with the case of prolate ellipsoids 
(£ = 1, dashed line) and of oblate ellipsoids (77 = 1, dashed-dotted 
line). In the inset the density is plotted as a function of 77 = £ for 
different polydispersities S r . N ~ 3000 — 4500 depending on the 
amount of the size polydispersity. 
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FIG. 6: The average coordination number as a function of the shape 
parameters r\ and £ . Here iV ~ 3000. 



Refs. |fl3l [l4ll where a slightly higher density is observed, 
due to their method for generating the packings: instead of 
being generated via deposition, their samples are prepared 
in order to be in jammed state using a generalized form of 
Lubachevsky-Stillinger algorithm ||28ll29l . 

The effect of polydispersity is examined by comparing the 
result for three different values of 5 r as shown in the inset 
of Fig. |5j>. The results indicate that the density seems to be 
insensitive to low and moderate size polydispersities. 

Figure [6] shows the average coordination number Z as a 
function of the shape of the particles. For spheroidal particles, 
Z increases gradually with the aspect ratios toward a maxi- 
mum. For a general ellipsoid the coordination number attains 
a maximum higher than the maximum observed for spheroids 
(77 = 1 or £ = 1), but decreases to values comparable to those 
of spheroids for larger 77 and £. This can be understood con- 
sidering the fact that for r\ = £ ^> 1, one of the semi-axis 
radii is significantly larger than the other two and therefore 
the particle can be regarded as a prolate ellipsoid. 



immediate realization is that the density is to a large extend 
symmetric with respect to the axis 77 = £. 

When the shape deviates from the spherical shape, i.e. when 
77 > 1 or £ > 1, the density first increases rapidly, then it at- 
tains a maximum and finally decreases slowly. From Fig. [5^ 
one sees that typically, along lines of constant 77 (or constant 
£) the maximum is observed around 77 = £. Further, along 
77 = £ a global maximum is attained at r\ = £ s» 1.4, as 
can be observed from Fig.[5]3 (solid line). This absolute max- 
imum curiously approximates the value of the golden ratio 
(1 + y5)/2, which from the definition of 77 and £ corresponds 
to when a = b + c. Whether such an observation is merely a 
coincidence or not needs further investigations which is out of 
the scope of our manuscript. 

These results are qualitatively similar to those obtained in 



B. Anisotropy 

To study the anisotropy of the deposits we investigate the 
orientational order of the contact vectors and the principal di- 
rections of the ellipsoids by calculating their angular distri- 
butions. In general two angles are needed to specify the ori- 
entation of a vector in three dimensional space, namely one 
azimuthal angle and one angle with the vertical direction. Our 
deposits are symmetric by construction in azimuthal direction. 
This is also observed in our results (not shown here for the 
sake of brevity). Therefore, we will only consider the angle 
with the positive vertical direction. 

Figure [7^ shows the results for deposits for different values 
of the shape parameters 77, £ where S r = . The vertical 
axis, P(fi(t9))df2, represents the fraction of number of vectors 
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FIG. 7: The anisotropy of the contact normal, (a) Histogram of the 
angle 9 between the contact vectors and the positive z-semiaxis for 
different shape parameters 77 and £, for the mono-sized case (S r = 
0). (b) Histogram of the angle 9 for for a given shapen (rj = 1.6, 
C = 1.8) but different size polidispersity. Here N ~ 3000 - 4500 
depending on the amount of the size polydispersity. 



within the solid angle swept by a conic shell whose axis is 
parallel to z and its vertex angle is 28. Note that the curves are 
symmetric with respect to 8 = ir/2 because for each contact 
normal of a particle there is a contact normal with opposite 
direction from the other particle to which it is in contact. 

The results show that for deposits of spheres the distribution 
P(n(8)) is independent of 9, and consequently the configura- 
tion of the contacts is isotropic and has no preferred direction. 
However, as the shape of the particles deviates from a sphere 
the number of side contacts between particles, i.e. the number 
of contacts on the horizontal plane, decreases and the contact 
normals line up along the vertical direction. This is seen in 
Fig. [7^ which shows for each shape there is a minimum at 
8 = 7r/2 which corresponds to the number of contacts on the 
(x, y)-plane, while having a maximum around 8 = (vertical 
direction). The minimum at 9 = tt/2 together with the maxi- 
mum at 8 = get more pronounced when increasing r\ and £. 
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FIG. 8: Distributions of the orientation angle 9 between the largest 
semi-axis a of each particle and the positive z-semi-axis. Each 
curves corresponds to particular values of the shape parameters. Here 
N ~ 3000. 



This means that for larger values of r\ and £ which correspond 
to higher asphericities the particles tend to lie horizontally, for 
minimizing potential energy during deposition. Consequently 
most of the contacts will occur with the particles beneath and 
above. This can have dramatic effect on the force chains and 
the response of under shear and compression. Such a study is, 
however, beyond the scope the current work. 

Figure|7^ also shows two crossing points at a value slightly 
higher than 8 ~ 7r/5, indicating this direction as an invariant 
directions while varying the shape of the depositing particles. 

Figure [T]} shows the comparison of the angular distribution 
of contact normals for deposits with 77 = 6,£ = 1.8, with 
different polydispersities 5 r . It is clear that the polydispersity 
has a minor effect. 

To study the anisotropy of the orientation of the particles 
we do a similar analysis for the principal directions of each 
ellipsoid. Figure [8] shows the results for the largest semi-axis 
a for the deposits. One can see that as the shape of the parti- 
cles deviates from a sphere the largest semi-axis a also tends 
to lie on the (x, y)-plane. An invariant orientation is also ob- 
served, now for 8 ~ tt/3. The polydispersities in the sizes of 
the particles has again no significant effects on the orientation 
anisotropy (not shown). 



V. CONCLUSIONS 

In this work, we studied systematically the static properties 
of deposits composed of ellipsoidal particles of different size 
and shape, focusing on their density and anisotropy. To this 
end, we developed an efficient Molecular Dynamics method 
for simulation of ellipsoidal particles which we used to gener- 
ate deposits of such particles under gravity. 

For mono-sized particles, the density increases rapidly as 
the shape of the particles deviates from a sphere, reaching for 
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particular values of the shape parameters a maximum signifi- 
cantly higher than the one observed in random close packing 
of spheres. With increasing further the shape parameters, that 
is, deviating further for a sphere density to decrease slowly. 
The densities observed in the deposits are generally lower 
than the ones obtained through procedures for optimal pack- 
ings. Thus, as expected, deposition is not an ideal process for 
achieving the largest densities. Typically, porous media and 
media formed through deposition processes present also den- 
sities much lower than the optimal one. Such observations 
point toward the possibility for using deposition of ellipsoids 
as a proper procedure for modeling porous media. This is the 
subject of future works. 

In the case of polydisperse particles the density takes higher 
values with qualitatively similar behavior. A rather unex- 
pected result is that the polydispersity in the sizes of the par- 



ticles has almost no effect on the anisotropic behavior of the 
system. 

We also study the anisotropy of such deposits and show that 
as the shape of the particles deviates from a sphere stronger 
anisotropics are observed. Such anisotropics appear due to 
the gravity and are absent in the systems without gravity. 
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